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Abstract: The application of the Markov chain to modeling agricultural succession is well known. 
In most cases, the main problem is the inference of the model, i.e. the estimation of the transi- 
tion matrix. In this work we present methods to estimate the transition matrix from historical 
observations. In addition to the estimator of maximum likelihood (MLE), we also consider the 
Bayes estimator associated with the Jeffreys prior. This Bayes estimator will be approximated by 
a Markov chain Monte Carlo (MCMC) method. We also propose a method based on the sojourn 
time to test the adequation of Markov chain model to the dataset. 
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Un modele markovien de dynamique d'usage des terres 

Resume : Les chaines de Markov sont depuis longtemps utilisees en modelisation de la dy- 
namique d'usage des terres. Dans la plupart des cas, se pose le probleme de I'inference du modele, 
c'est a dire de la construction de la matrice de transition qui dirige la dynamique de succession. 
Nous presentons dans cet article des methodes pour estimer cette matrice a partir d'un historique 
d'observations. En plus de I'estimateur du maximum de vraisemblance (EMV), nous considerons 
I'estimateur bayesien associe a la loi a priori non informative de Jeffreys. Cet estimateur bayesien 
sera approche par une methode de Monte Carlo par chaine de Markov (MCMC). Nous etudions 
egalement I'adequation entre les temps de sejour, en un etat, constates dans les donnees et leur 
estimation par le modele de Markov. 

Mots-cles : Modele de Markov, Monte Carlo par Chaine de Markov, loi a priori de Jeffreys, 
dynamique d'usage des terres. 
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1 Introduction 

Population pressure is one of the major causes of deforestation in tropical countries. In the region 
of Fianarantsoa (Madagascar), two national parks Ranoniafana and Andringitra are connected 
by a forest corridor, which is of critical importance to maintain the regional biodiversity. The 
need for cultivated land pushes people to encroach on the corridor to look for swallows to be 
converted into paddy fields, and then to clear slope forested parcels for cultivation. Once swallows 
are all converted in paddy fields, the dynamic of slash and burn cultivation is clearly opposed 
to the dynamic of forest conservation and regeneration. To reconciling forest conservation with 
agricultural production, it is important to understand and model the dynamic of post-forest land 
use of these parcels. We will use a first data set developed by IRD, in the western edge of the 
corridor, consisting of the annual state of 43 parcels initially in forest, during 22 years since first 
clearing (Figure 1). Each parcel can take four possible states: forest (F), annual crop (C), fallow 
(J), perennial crop (B). 

The use of Markovian approaches to model land-use transitions and vegetation successions is 
widespread [11, 12, 13]. The success of these approaches is explained by the fact that agro-ecological 
dynamics are often represented as discrete succession of a finite number of states, each one with 
its holding time. Both agronomists and ecologists, in dialog, actually fail in predicting the future 
succession of these states, knowing the previous land-use history. They ask the mathematicians for 
detailing the characteristics of these dynamics and defining how to pilot them. The construction, 
manipulation and simulation of such models are fairly easy. The transition probabilities of the 
Markovian model are estimated from observed data. The classical reference [1] proposes the maxi- 
mum likelihood method to estimate the transition probabilities of a Markov chain. An alternative 
is to consider Bayesian estimators [10, 8[. In this paper we explore and test several modeling tools, 
Markov chain, Bayesian estimation and MCMC procedure to better fit with the actual data. These 
results are needed by agro-ecologists who try to model the land use dynamics, at a parcel scale. 

The model is introduced in Section 2, then the maximum likelihood estimator and the Bayesian 
estimator are presented in Sections 3 and 4 respectively. These estimators are applied to simulated 
data in Section 5 and to the real data set in Section 6. The Markov model is evaluated in Section 
7. Conclusion and perspectives are drawn in Section 8. 

2 The model 

We make the following hypothesis: 

(Hi) The dynamics of the parcels are independent and identical. 

This means that {en)n'^o-2i ^^^ ^^ independent realizations of a same process {Xn)n=0:2i- This 
assumption is not realistic as the dynamics of a given parcel depends on: 

• farmer decisions; 

• exposition, slope and distance from the forest, that means properties of the same plot; 

• neighboring parcels. 

This assumption, however will lead to a simple model. 
We also suppose that: 

{H2) The process {Xn)n=o---2i is Markovian and time-homogeneous. 

The homogeneity assumption is also simplistic but we assume that the transition law of parcels 
will poorly varied during this 22 year period. 
Finally we suppose that: 
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parcel nuniber 
12 3 4 ... 



Figure 1: Annual states (e5^)^^o-2i corresponding to ^3 parcels and 22 years. These parcels are 
located on the slopes and lowlands on the edge of the forest corridor of Ranomafana-Andringitra, 
Madagascar. The states are: forest (F), fallow (J), annual crop (C), perennial crop (B). 




B 1 

perennial crop 



Figure 2: Four states Markov chain diagram: forest (F), annual crop (C), fallow (J), perennial 
crop (B); F is the initial state and B is an absorbing state. 



{H^) The initial state is F . 

These hypotheses lead to a model X = {Xlf=J-Q^^_^, P = 43, N = 22 where (XP)„^o...Ar_i 
are P independent Markov chains, with initial law dp and transition matrix Q of size 4x4. The 
state space is: 

E'^{F,C,J,B}. 
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Hence: 

p 

p=l 

= f[5F (eg) Q(eg, e?) • • • Qie^^_„ e^^__,) (la) 

for all eP, £ £;. 

Some transitions are not observed at least in the time scale considered here: once the parcel 
leaves the state "forest" by first clearing, it cannot come back; and similarly when it reaches the 
state "perennial crop", it stays there during the sample time, that means permanently in this model. 

{H4) The transitions C^F,J^fF,B^^F,B^^C,B^^J,J^^Bdo not exist in the model, 
all other transitions are possible. 

In particular: once the parcel leaves the state "forest", it cannot come back; when it reaches the 
state "perennial crop", it stays there permanently. 

This hypothesis implies that (i) the realistic transitions, which are not observed during the 
considered time scale, J ^- F, i? — >■ C, do not exist in the model, {ii) the unrealistic transitions 
C—>-F,B—>-F,B—>-J and J —>^ B do not exist, {Hi) the state F is transient (more precisely when 
the chain leaves the state F it will never come back to that state), (iv) the state B is absorbing. 

To summarize we consider a transition matrix Q of the form: 

[1-01-02 Oi 02 \ 

1 — P3 — P4 ^3 P4 

05 I-6I5 

\ 1/ 



Q = 



(lb) 



that corresponds to Figure 2, and depends on a 5-dimensional parameter: 

= (^1,^2,^3,^4,^5) 
belonging to the set: 

e = {0 e [0, 1]5 ; 0i + 02<l,03 + O4<l}- (2) 

Let Pg denotes the probability under which the Markov chain admit Q with parameter as a 
transition matrix. 

3 Maximum likelihood estimation 

We recall the classical results of Anderson-Goodman [1] to compute the MLE of the matrix Q. 
The likelihood function associated with {Pg; £ Q} is: 

Li0) ^='- PgiX'oiN-i = el^:N-i) = n 5f (eg) g(eg, e{) ■ ■ ■ Qie%^,, e^_,) 

p=i 

where Q is defined by (lb), for any eQ:^_j^ G E-^^^ . Let n^^, be the number of transitions from 
state e to state e' for a parcel p in X: 

"L' = XI l{^S-i=e} l{xs=e'} ye,e'eE (3) 
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and Tiee' be the total of number of transitions from state e to state e' (e, e' Cz E): 



p 

dot 



E<e'- (4) 

p=l 



According to (lb): 



L{e) = Yl Q(l, 1)""^^ Q(l, 2)<c Q(i, 3)""^./ Q(2, 2)"cc Q(2, 3)"&^ 

Q(2,4)"&B Q(3,2)"5cQ(3,3)"5. 
p 
= n^^ " ^1 ~ ^2)"^^ 6'^^"' C^' (1 - 6*3 - 6'4)"&c 6»3&-^ 6'4^« (1 - e^)""-'-' el^"' 
p=l 

and from (4): 

L{e) '^ (1 - 6I1 - 6(2)"^^ 6'5'^c 6i^^-^ (1 - 6*3 - 6'4)"'='= 0^^'^ e^""" (1 - ^s)"'-' C'^ 

so that the log-likelihood function reads: 

l{e) = npF log(l - 6*1 - 6*2) + n^^c log(6'i) + npj \og{92) 

+ ncc log(l - 6*3 - 6*4) + ncj log(6'3) + ncB log(6'4) 

+ njjlog(l-6l5)+njclog(6'5). (5) 

The MLE 9 is solution of dl{9)/d9\g^g = 0, that is: 

dlje) 
d9i ~~ 
dl{e) 



903 

dl{6) 
We get: 



npF npc 


= 0, 
= 0, 


dl{6) 
392 

dl{9) 
894 


npF np.i 


{1-61-62) 61 
ncc ncj 


[1-61-92) 92 
ncc ncB 


(1-6,-64) 6, 
njj njc 


{1-6,-94) 64 


(1-6,) 6, ''■ 





= 0, 
= 0, 



npC n ^FJ 



npF + npc + npj npp + npc + npj 

ncj a "^CB 
, t'4 = 

ncc + ncF + ncj ncc + ncj + ncB 
njc 



75 



njc + njj 

4 Bayesian estimation 

We suppose that an a priori distribution law 7r{6) on the parameter 9 is given. According to the 
Bayes rule, the a posteriori distribution law 7r{9) on 6 given the observations X is: 

^(61) ex 1(6) Tr{6) (6) 

where L{6) is the likelihood function. The Bayes estimator 9 of the parameter 6 is the mean of 
the a posteriori distribution: 

6^^=^'' f en{9)d6. (7) 

Je 
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4.1 Jeffreys prior 

Numerical tests that will be performed in Section 5.1 suggest that the Jeffreys prior is well adapted 
to the present situation and we introduce it now. This prior distribution (non-informative) is 
defined by [7]: 

7r(6i) oc y/det[I{6)] 



where I{9) is the Fisher information matrix given by: 



E, 



1(9) = 
and l{d) is the log-likelihood function. Hence 

1(0) 

with 



V de.ddJ 



l<k,l<5 



Ai.2 

^3,4 

a5 



Ak,e = -Eg 






\ dBkdei 

So deil{9) = det Ai,2 x det ^3^4 x 05 and 



an(9) 
ri{ 






05 



'Eff 



d^l{9) 



t:{9) cx ^detAi.2 x detyl3_4 x 05. 



According to (5): 



3^;(e 



dH(e_ 

d^92 

d^l(9 



d^9:: 



d^9i 

d^l(9 
d^95 



{1-61-62)^ '^ 9j 
(l-ei-82)^ '^ 91 



{l-63-eiY ^ 9l 



(1-63-64)2 ^ ' 
rijj 1 njc 

(i_e5)2 -I- ei 



and 



detAi,2 = (pS 



"ff] 



Ee["Fc] 



-92) 



detA 



05 



lEs[rtjj] I Ee[njc] 
(1-05)2 ^ e| 



\{l-9i- 



_ f E(,[kcc1 I Eelncj] \ ( ^ 

3,4 - (^(i_e3_e^)2 + 02 j (^(1 



2)2 



Ee[Ttcc] 
-63-64)2 



d^l{9) 



a6i962 " 


(1- 


-61-62)2 ' 


a2;(6l) 




npF 


902 96i 


(1- 


-61-62)2 ' 


a2;(0) 




ncc 


963 964 


(1- 


-63-64)2 ' 


d^m 




ncc 


904 963 


(1- 


-63-64)2 ' 



EeKvl 



Eejncs] 



V(i-ei-02)2y 

V (1-63-64)2 ; 



(8) 



(9a) 



(9b) 

(9c) 
(9d) 
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From (3) and (4): 

P N 



EeKe'] = ^^Pe(XP = e',XP_, = e) 

p=l n=2 

N 

= P ^Pe(X„ = e'|X„_i - e)P0(X„_i = e) 

11=2 n=2 

AT 

= PQ(e,e')5][Q("-^)](P,e) (9e) 



for all e,e' e £'. 

Note that [2] proposed a more complex method to compute the Jeffreys prior distribution. 

4.2 MCMC method 

Although the Jeffrey prior distribution is explicit, we cannot compute analytically the correspond- 
ing Bayes estimator. We propose to use a Monte Carlo Markov chain (MCMC) method, namely a 
Metropolis-Hastings algorithm with a Gaussian proposal dsitribution, see Algorithm 1. 

choose 9 

for fc = 2, 3, 4, . . . do 

e~AA(0,a2) 

u - C/[0, 1] 

if u < a then 

9 <- 9""°" % acceptation 
end if 



k k 

end for 



Algorithm 1: MCMC method: Metropolis-Hastings algorithm with a Gaussian proposal distribution. 
The target distribution is 7t{9) defined by (6), the Gaussian proposal PDF (probability density 
function) is g(- — 9) (g PDF of the A/'(0,cr^) distribution) where 9 is the current value of the 
parameter. 
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True values 

— - Bayesian estimates with Jetfreys prior 
• — I\^aximum iil<eliliood estimates 

■ ■ ■ Bayesian estimates with beta prior 

— Bayesian estimates with uniform prior 





Estimation of p 



Estimation of q 



Figure 3: Two states case with paramaters p = P(X„+i = 0|X„ = 0) and q = V{Xn+i = l|Xn = I)-' 
we compare the following priors: uniform, beta(l/2, l/2j and Jeffreys. Jeffreys prior gives better 
results than the two other priors. 

5 Simulation tests 

5.1 Two states case 

We first consider the simpler two states case E = {0, 1}. It has no connection with the Markov 
model considered in the present work but it allows to easily compare the following different prior 
distributions: 

(i) the uniform distribution; 

(ii) the beta distribution of parameter (^ i); 

(Hi) the non-informative Jeffreys distribution. 

The Bayesian estimator is explicit for the two first priors, see Appendix A. 

We compare the MLE and the Bayesian estimator with uniform and beta priors, that can 
be explicitly computed, see Appendix A, with the Bayesian estimator with Jeffreys prior that is 
computed by an MCMC method that will be explained later. Results proposed in Figure 3 tend 
to demonstrate that the Jeffreys prior gives better results than the two other priors. 

5.2 Four states case 

Before processing the real data set of Figure 1 with a four states Markov model, we consider a 
simulated case test. We aim to compare the MLE and the Bayes estimator with the Jeffreys prior. 
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Bayesian estimates 
Maximum iil^elihood estimates 
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Figure 4: Empirical PDF for the error terms (12) associated with the 2-norm (10) and the Frobenius 
norm (11) based on 1000 simulation of the parameter 9. 

We compute the distance between the real transition matrix Q and its estimation, with the 
MLE or the Bayes estimator, given by the Frobenius norm: 



|^||2,^=trace(A*yl) 



and the 2-norm: 



||a||2=Va._(a*a) 



(10) 



(11) 



where \^^^{A* A) is the largest eigenvalue of the matrix A* A. 

We sample 1000 independent values 6*'^^^ of the parameter according to a uniform distribution 
on Q defined by (2), that is a uniform distribution on [0, 1]^ with specific the constraints. For each 
H., we simulated data (eJ^)fjZQ'.2i according to the model (1) that is with the transition matrix Qgt 
defined by (la). Then we compute the MLE ^'^^ and the Bayes estimate ^'^^■' with the Jeffreys 
prior. Then we compute the errors: 






(12a) 
(12b) 



for the two different norms. 

In Figure 4 we plotted the empirical distribution of the errors e^ and e^ , d = 1- ■ ■ 1000, for the 
two different norms. We see that the Bayes estimator give slightly better results than the MLE. 

6 Application to the real data set 

For the real data, the result of both approaches are slightly different (see Figure 5 and the Table 
1). This calls into question the considered model. We develop this point in the next section. 
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Bayesian estimates 

— Maximum lil<elihood estimates 










D.QO D.01 Q.02 0.03 0.04 0.05 




Figure 5: MLEs and Bayes estimators for the real data set: posterior empirical distributions given 

by the MCMC iterations ( — /green) and the associated mean ( /red) and the maximum likelihood 

estimates (■ ■ ■ /blue). 
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Bayesian estimates 
to 





F 


C 


J 


B 1 


F 


0.9121 

IT 


0.0842 
0.7417 


0.0037 


^H 


C 


0.2433 


0.0150 


J 


0.3273 


0.6727 


■■ 


B 






^^ 



Maximum likelihood estimates 
to 





F 


C 


J 


B 1 


F 


0.9158 


0.0823 


0.0019 ^^^H 


C 


[:■ 


0.7449 


0.2426 


0.0125 1 


J 


0.3233 


0.6767 


^^^H 


B 







1 



Table 1: Bayesian and Maximum likelihood estimates. Dark grey cells correspond to transition 
probabilities that are supposed to be known; light grey cells correspond to transition probabilities 
that are deduced from the other ones. 

Distribution of the time to reach B 

Given the two estimations of the transition matrix, we would like to address the two following 
questions. First, what is the distribution law of the first time to reach the absorbing state B ? 
Second; as B is absorbing and all other states are transient, the limit distribution of the Markov 
chain is 5b , but before this state b is reached what is the "limit" distribution of X„ on the other 
states ? This distribution is called the quasi-stationary distribution of the process X„ and we will 
compute it. 

To answer the first question we use the result of Appendix B : the distribution law of the first 
time TpB to reach B starting from F: 



{tfb = n\Xa = F)= P(X„ = B, Xr,, ^ B, m = 1, 



l\Xo = F) 



is given by recurrence formula (16) and plotted in Figure 6 for both the Bayesian and the maximum 
likelihood estimates. The mean time is 92 years for the Bayesian estimate and 96 years for the 
MLE. 

Limit distribution before reaching B (quasi- stationary distribution) 

The answer to the second question is given by the so-called quasi-stationary distribution, see 
Appendix C. From (17) we can compute the quasi-stationary distribution ft — {fl{F) , il{C) , fi-{J)) 
associated with the estimators of Q with the maximum likelihood and the Bayesian approaches. 
The results are: 

Bayesian estimator 
Maximum likelihood estimator 

Hence, conditionally the fact that the process does not reach B, and as soon at it leaves the state 
F, it will spend 57% of its time in the C state and 43% of its time in the J state. 



F 


C 


J 





0.5659 


0.4341 





0.5672 


0.4328 
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Bayesian estimate 

— - Maximum lil<eliliood estimate 




Figure 6: Distribution of the time to reach the state "B". The mean time is 92 years for the Bayesian 
estimate and 96 years for the MLE. 
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7 Model evaluation 

In this section we test the fit between the data and the model. From the data set of Figure 1 it is 
clear that the holding time in the state "forest" does not seem to correspond to that of a Markov 
chain. The holding time S{e) of a given state e G E, also called its sojourn time, is the number of 
consecutive time periods the Markov chain Xn remains in this state: 

S{e) "^ inf{n e N; X„ ^ e} 

conditionally on Xq — e. The distribution law of S{e) is given by: 

Pe(5(e) = n) = Pe(Xo = • • • = X„ = e, X„+i ^ e) 

= Y, ^<^(^o = ■■■ = X„ - e, X„+i = e') 

= J2 Q(e,e)- - -Q(e,e) Q(e, e') = {Q{e, e))" (1 - Q(e, e)) 

for n > 1 and for n = 0, that is a geometric distribution of parameter Q{e,e) — V{Xn+i = 
e\Xn = e). Note that Ee5'(e) = 1/(1 - Q{e,e)) and vare(i?e) = <9(e,e)/(l - Q{e,e)f. 

7.1 Goodness-of-fit test 

In order to test if the distribution of the holding time S{e) on each state e G E of the data set 
(6n)^=o-2i i^ geometric, we use a bootstrap technique for goodness-of-fit on empirical distribution 
function proposed in [5]. 

Considering a sample Si, . . . ^ Sk of size k from a discrete cumulative distribution function 
{F{n))neN, we aim to test the following hypothesis: 

Ho : Fe¥ = {Fp-.peS}. (13) 

In our case, Fp is a geometric cumulative distribution function (CDF) with parameter p G [0 1]. 
Classically, we consider an estimator: 

of p and we compute the distance between the theoretical CDF Fp and the empirical CDF: 

d f 1 '^ 

We use the Kolmogorov-Smirnov distance defined by: 

K* = K{Si.,k) = sup Vk\Fs,Jn) - i^T(5,-.)(«)l (14) 

To establish whether K* is significantly different from or not, we simulate M samples of size k: 



and we let: 



dot 



where K is the function defined in (14). 
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F (forest) 


Holding time values 
Number of occurrences 


1 
1 


3 
9 


11 
2 


13 
1 


14 
6 


15 
21 


16 
3 


C (annual crop) 


Holding time values 
Number of occurrences 


1 
11 


2 

17 


3 

12 


4 
5 


5 
9 


6 

7 




J (fallow) 


Holding time values 
Number of occurrences 


1 
16 


2 
12 


3 

7 


4 
4 


6 
2 


8 
1 


11 
1 



Table 2: Holding time values (year) on each state F, C, J in the data set. 



The p- value associated to that test is: 



M 



dcf 1 v-^ 



If p is less than a given threshold a, corresponding to the probability chance of rejecting the null 
hypothesis Hq when it is true a tort, then Hq is rejected. 

7.2 Holding time goodness-of-fit test 

The B state (perennial crop) is absorbing so that its holding time is infinite. Moreover, states that 
appear at the end of the series of Figure 1 are not treated (they are considered as censored data). 
Then the holding time values on each state F, C, J in the data set are given in Table 2. 

In order to test the hypothesis H^ we use the MLE for the parameter p of the geometric PDF: 

(15) 



1 + i E£=i St 



Indeed, the likelihood function is: 



L(p) = (l-p)^>---(l-p)^'=p 

and L'{p) = leads io n — p (X^fci S^ + n) = Q and (15). 
The complete test procedure is given by Algorithm 2. 



(i_p)5:tiS.p'= 



7.3 Results 

In Figure 7 we plotted the empirical PDFs of holding time of states "Forest", "Annual crop" and 
"Fallow", associated with the data set of Figure 1, and the geometric PDF corresponding to the 
parameter p estimated by (15). We see that in the case of the "Forest" state the matching is 
questionable. In Figure 8 we plotted the empirical PDF of K™ for the three states. The value of 
K* and of the associated p-values are: 





"Forest" F 


"Annual crop" C 


"Fallow" J 


K* 


3.051 


1.086 


1.104 


p- value 





0.224 


0.255 



In conclusion, the geometric distribution hypothesis is strongly rejected for the state "forest". 
The p-value for this state is null. This is understandable as this state does not really correspond 
to a "dynamic state". 
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n ^ sup(S'i,.^. ,5*^;) 

K* <~ svLp{Vk\F'^^'-''{n) - i^^('5i^'')(n)| , < n < n} 

for TO = 1,2, . . . ,M do 

cm Q771 ^^^ p 

Di , . . . ,Di^ ^p 

K"" ^ supiVklF^^i-in) - F'^(^"'=)(n)| , < n < n} 

1 v-^A/ -, 
P^ Jl l^rn=l ^{K'">-R'-} 

end for 

if p < a then 

accept Hq, 
else 

reject Hq. 
end if 



Algorithm 2: Parametric bootstrap for goodness- of -fit with the geometric distribution of parameter 
P- 





10 20 30 40 50 60 

Annual crop C 




10 20 30 

Fallow J 



Figure 7: Empirical PDFs of holding time of states {F, C, J) associated with the data set of Figure 
1; and the geometric PDF, represented as continuous read lines, corresponding to the parameter p 
estimated by (15). 




1 2 

Forest F 




0.5 1.0 1.5 2.0 

Annual crop C 



J 



I r r 



0.5 1.0 1.5 2.0 2.5 

Fallow J 



Figure 8: Empirical PDF of K"^ (sample from K* ) for the three states {F, C, J) et K* value 
(vertical line). 
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8 Conclusion and perspectives 

We proposed a Markovian model of land use dynamics for parcels near the forest corridor of 
Ranomafana and Andringitra national parks in Madagascar. We supposed first the dynamic uses 
of the parcels are independent and identically distributed; second that the dynamics is Markovian 
with four states. The transition matrix depends on five unknown parameters. We considered the 
MLE and the Bayes estimate with Jeffrey prior. In this last case, the estimator is computed with a 
MCMC procedure. The Bayes estimator performs slightly better than the MLE. On the real data 
set, the two estimators give rather similar results. 

We assessed the adequacy of the model to real data. We focused on the holding times: we 
tested if the empirical holding times correspond to a geometric distribution. We used a parametric 
bootstrap goodness-of-fit on empirical distribution. Clearly the geometric distribution hypothesis 
is violated in the case of the "Forest" state. 

The "Forest" state therefore requires a special treatment. In a near future we are now developing 
a semi-Markov model where the sojourn time on the state F will better match the data set and so 
will not be geometric. 

The long time behavior of the inferred model is dubious as the present data set is relatively 
limited in time (22 years). This data set implies a relatively short time scale where some rare 
transitions, like the forest regeneration, are not observed. Note that the Bayesian approach has an 
advantage over the likelihood approach in that it allows to incorporate prior knowledge about these 
rare and unobserved transitions. The likelihood approach will set their probabilities zero while the 
Bayesian approach will incorporate a priori knowledge and assign them positives probabilities. A 
new database is currently being developed by the IRD. It will be for a longer period of time and 
a greater number of parcels, it will also allow to consider a more detailed state space comprising 
more than four states. In a longer time scale, it is reasonable to suppose that F and B have long 
sojourn time distributions, the one associated to F being longer than the one associated to B. 
Also B will not be absorbing anymore as well as the forest regeneration will be possible, i.e. the 
transition from J to F will be possible. The associated model will present multi-scale properties, 
namely slow and fast components in the dynamics, that will be of interest. 

Part of the complexity of these agro-ecological temporal data comes from the fact that some 
transitions are "natural" while others come from human decisions (annual cropping, crop abandon- 
ment, planting perennial crops, etc.). It should also be interesting to study the dynamics of parcels 
conditionally on the dynamics of the neighbor parcels. This model could be more realistic but re- 
quires first studying the farmers' practices in order to limit the number of unknown parameters in 
the model. 
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Appendices 

A. Explicit Bayes estimators for the two state case 

Let {Xn)Q<n<N be a Markov chain with two states {0, 1} and transition matrix 

We suppose that the initial law is the invariant distribution fi = {q/p + q,p/p + q), that is the 
solution of /iQ = /i. The unknown parameter is 9 = {p,q) G [0, 1]^ and the associated likelihood 
function is 

Ln{0) = VeiXo-.N = xo:n) = (1 - p)""" p""^ g"^" (1 - g)"" 

where riij = nij{xQ-N) — X]n=o ^{x„=i} l{x„+i=j} is the number of transition i — > j in xo:jv- 

We consider the following priori distributions: the uniform distribution tt" on [0,1]^ and the 
beta distribution tt^ with parameters (a, &), that is 

1 



p{a,b) 



where /3(a, 6) is the beta function: 

Jo r(a + b) 

with r(z) = /g+°° i^-i e"* dt. Here we will choose a = b= 1/2, note that r(l) = 1 and r(i) = ^tt. 
For these two priors we can explicitly compute the posterior distriution and the associated Bayes 
estimators. Indeed the posterior distribution n is given by the Bayes formula: Tf{9) ex Lf^{9) Tr{9), 
that is: 

if"(6l) (X Lpfie) = (1 - p)"™ p""i q"i" (1 - g)"" , 

^^(61) ex LN{e) 7r^(6l) = (1 -p)"™-5 /'Ol-i ^"10-1 (1 _ g)«ll-5 

and the corresponding Bayes estimator are: 



en''{e)d9, 61^ == / 9n''{9)d9. 

[0,1]^ J[oa]2 

We can easily check that the estimators of p and q for the uniform prior: 

ngi + 1 
"-01 + "00 + 2 ' 



-U noi +1 .u "10 + 1 

P = ] tt: ' 1 ~ 



and for the beta prior: 

p-= "°^ + ^ , ,- 

"01 + "00 + 1 

Note that in this case the MLE estimators are: 

^MLE "Ol ~EMV 

"00 +"01 ' "11 +"10 
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B. Distribution law of the time to reach a given state 

Let Xn be an homogeneous Markov chain with finite state space E and transition matrix Q. We 
aim to get an explicit expression of the distribution law f^", = V{Tee' = n), n e N, of the first 
time Tee' to reach state e' after leaving state e defined as: 

Tee' = inf{n > 1 : X^^ e'\Xo = e}. 

For n > 1: 

Q(")(e,e')=P(X„ = e'|Xo = e) 

= P(X„ = e', Tee' = 1|^0 = e) + P(X„ = e', Tee' = 1|^0 = c) + • • • 

• • • + P(X„ = e', Tee' = n - l|Xo - e) + P(X„ = e', r^e' = n\Xo = e) 
= /iVg^-'He', e') + /(f) Q("-2)(e', e') + • • • + /,'r '^Q^'^(e', e') + /i"? 

hence /^gV could be computed recursively according to 

n-l 

/i") =Q(")(e,e')-5]/i'?Q("-'=He',e') (16) 



fc=i 



with/iV=Q(i)(e,e'). 

C. Quasi-stationary distribution 

We consider the probability to be in e G {F, C, J} before reaching i? and starting from F: 

finie) = P(X„ = e|X,„ ^ B, m = 1, . . . , n - 1, Xo = F) 
= P(X„ = e|X„_i ^ S, Xo = F) 
^ PjXn = e, Xn-i ^ -B|Xo ^ f ) 
P(X„_i ^ S|Xo = F) 
P(X„ = e|Xo = F) 



When 



^-j:e'e{F,c,J}nXn-i=e\Xo = F)- 



lin{e) > A(e) , e e {F, C, J} 



the probability distribution {iJ-{e))e£{F,C,J} is called quasi-stationary probability distribution. This 
problem was originally solved in [3]: fi = [fi{F) JJl{C) p-{J)] exists and it is given by the equation 

fiQ^~Xfi (17) 

with jj,{e) > and /i(F) + fl{C) + jl{J) — 1, where Q is the 3x3 submatrix defined by 



Q=' ^ 





and A is the spectral radius of Q. 
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